Customized variables

The analysis depends on the following variables defined in configure.R

Variable Value
Experimental conditions rbohC, rbohF, WT
Data file(s) data/rbohC_C_vs_Cd.tsv, data/rbohF_C_vs_Cd.tsv, data/WT_C_vs_Cd.tsv
Total proteins file data/identification_all_samples_proteins.tsv
Subcellular location file with annotations location/proteins_locations_modified.tsv
Subcellular location ontology location/subcel_location_unique.txt

Load data files

Read proteomic results of DEPs for different analyses.

# DEP counts ####
counts <- list()
for (i in 1:length(FILE)) {
  counts[[FILE_NAME[i]]] <-  read.table(file = FILE[i],
                                        sep = "\t",
                                        header = TRUE,
                                        row.names = 1,
                                        check.names = FALSE,
                                        comment.char = "") # turn off the interpretation of comments
}

# all proteins ####
all_proteins <- read.delim(file = FILE_TOTAL,
                           sep = "\t",
                           header = TRUE,
                           row.names = 1)

# subcellular location file ####
location <- read.delim(file = LOCATION,
                       sep = "\t",
                       header = TRUE,
                       row.names = 1)

# subcellular location ontology ####
location_cat <- scan(file = LOCATION_CAT,
                     what = character(),
                     sep = "\n")

# construct final datatable
final_datatable <- data.frame()
# merge data from counts
for (i in 1:length(counts)) {
  tmp <- counts[[i]][, grep("Abundance Ratio \\(log2\\)", colnames(counts[[i]])), drop = FALSE]
  colnames(tmp)[1] <- paste0("log2_", names(counts[i]))
  tmp$Protein <- rownames(tmp)
  tmp <- tmp[, c(2,1)]
  if (length(final_datatable) == 0) {
    final_datatable <- tmp
  } else {
    final_datatable <- merge(final_datatable,
                             tmp,
                             by = 1,
                             all = TRUE,
                             sort = FALSE)
  }
}
# merge data from location
final_datatable <- merge(final_datatable,
                         location,
                         by.x = 1,
                         by.y = "row.names",
                         sort = FALSE)
final_datatable <- final_datatable[, c(1, 7, 8, 2, 3, 4, 5)]
final_datatable$Transmembrane <- gsub("TRANSMEM.*", "TRUE",
                                      final_datatable$Transmembrane)
final_datatable[is.na(final_datatable)] <- 0

Heatmap

Expression of DEPs in WT, rbohC and rbohF as normalised abundances.

heatmap_list <- list()
col_min_max <- c("white", "#4292C6")
my_palette <- colorRampPalette(col_min_max)(50)
for (i in 1:length(counts)) {
  tmp <- counts[[i]]
  # keep only abundances columns
  tmp <- tmp[-c(1, 8, 9)]
  # change colnames
  cols <- colnames(tmp)
  cols <- gsub(".*Sample..", "", cols)
  colnames(tmp) <- cols
  name <- names(counts[i])
  tmp <- as.matrix(tmp)
  heatmap_list[[name]] <- tmp
  distance <- dist(tmp, method = "euclidean")
  tmp_hc <- hclust(distance, method = "ward.D2")
  heatmap(tmp,
          Colv = NA,
          Rowv = as.dendrogram(tmp_hc),
          scale = "row",
          cexRow = 0.5, # nombres pequeños para que quepan
          cexCol = 0.8,
          col = my_palette, # va de 0 a infinito, no hay negativos
          main = paste0(names(counts[i]), " (", nrow(counts[[i]]), " DEPs)")
  )
  legend(x = "topright", 
         inset = c(-0.01, 0),
         legend = c("Low", "High"),
         cex = 0.8, 
         fill = col_min_max,
         bty = "n"    # no box around
  )
}

Mutant comparison

We compare the differentially expressed proteins (DEPs) that are common between WT, rbohC and rbohF.

DEP intersections: Venn diagram

# list to save the ids of DEP
proteins_ids <- list()

# DEPs including up- and down-regulated (All-DEPs)
# The first level of the list are the rbohC, rbohF and WT genotypes
proteins_ids[["all"]] <- list()
for (i in 1:length(counts)) {
  tmp <- rownames(counts[[i]])                # get only names
  proteins_ids$all[[names(counts[i])]] <- tmp # assign protein ids to genotypes
}

# up- and down- DEPs split ####
proteins_ids[["up"]] <- list()   # empty list for UP
proteins_ids[["down"]] <- list() # empty list for DOWN

# recorrer los tres elementos de la lista (los tres genotipos)
# y dividir los DEP entre UP y DOWN
for (i in 1:length(counts)) {
  # Qudarse con el valor de la columna que dice si es una DEP UP o DOWN
  tmp <- counts[[i]][, grep("Abundance Ratio \\(log2\\)", colnames(counts[[i]])), drop = FALSE]
  tmp_up <- tmp[tmp[1] > 0, , drop = FALSE]       # seleccionar los UP. drop=FALSE para conservar nombre
  tmp_down <- tmp[tmp[1] < 0, , drop = FALSE]     # seleccionar los DOWN
  name <- names(counts[i])                        # nombre de genotipo para la lista
  proteins_ids$up[[name]] <- rownames(tmp_up)     # subllista de UP para cada genotipo
  proteins_ids$down[[name]] <- rownames(tmp_down) # sublista de DOWN para cada genotipo
}

# venn diagrams for all-DEPs, up and down ####
venn_intersections <- list()       # para guardar las intersecciones para tabla de discordancias
venn_data <- data.frame()
for (i in 1:length(proteins_ids)) {
  ## Venn #####
  venn_plot <- venn(proteins_ids[[i]],
                    ilabels = "counts",
                    zcolor = "style",
                    box = FALSE)
  title(paste0(names(proteins_ids[i]), "-DEPs"),  # toma el nombre del id de la lista
        line = -1)                                # para que no se pegue a la imagen
  if (length(venn_data) == 0) {                   # Añadir el primer elemento
    venn_data <- venn_plot[2:8, 4, drop = FALSE]
    venn_data$intersections <- rownames(venn_data)
    venn_data <- venn_data[, c(2,1)]
    colnames(venn_data)[i+1] <- names(proteins_ids[i])
  } else {                                        # merge los demás elementos por ID de prot
    tmp_venn <- venn_plot[2:8, 4, drop = FALSE]
    tmp_venn$intersections <- rownames(tmp_venn)
    tmp_venn <- tmp_venn[, c(2,1)]
    colnames(tmp_venn)[2] <- names(proteins_ids[i])
    venn_data <- merge(venn_data, tmp_venn, by = 1, sort = FALSE)
  }
  
  intersections <- calculate.overlap(proteins_ids[[i]])
  ## define names for the intersections ####
  a123 <- paste0(names(proteins_ids[[i]][1]),
                 "-",
                 names(proteins_ids[[i]][2]),
                 "-",
                 names(proteins_ids[[i]][3]))
  a12 <- paste0(names(proteins_ids[[i]][1]),
                 "-",
                 names(proteins_ids[[i]][2]))
  a13 <- paste0(names(proteins_ids[[i]][1]),
                 "-",
                 names(proteins_ids[[i]][3]))
  a23 <- paste0(names(proteins_ids[[i]][2]),
                 "-",
                 names(proteins_ids[[i]][3]))
  a1 <- paste0(names(proteins_ids[[i]][1]))
  a2 <- paste0(names(proteins_ids[[i]][2]))
  a3 <- paste0(names(proteins_ids[[i]][3]))
  
  names(intersections) <- c(a123, a12, a13, a23, a1, a2, a3) # give the new names to intersections
  # asingar los solapamientos entre fenotipos al elemento de la lista all, up o down
  venn_intersections[[names(proteins_ids[i])]] <- intersections
}

Save DEP intersections

filenames <- ""  # empty variable for filenames
for (i in 1:length(venn_intersections)) {
  ## modify intersection results ready to save
  tmp <- venn_intersections[[i]]
  tmp <- do.call(rbind, lapply(tmp, data.frame)) # convert list in table
  tmp$intersections <- rownames(tmp)                  # add a column with the comparisons
  colnames(tmp)[1] <- "protein"                 # change column name to 'protein'
  tmp$intersections <- gsub("\\..*", "", tmp$intersections) # remove the .Number
  tmp <- tmp[,c(2,1)]                           # change column order
  file_name <- paste0(SUB_DIR,
                      "DEP_",
                      names(venn_intersections[i]),
                      ".tsv")
  write.table(tmp,                              # save the table
              file = file_name,
              quote = FALSE,
              sep = "\t",
              row.names = FALSE,
              col.names = TRUE)
 filenames <- c(filenames, paste0("**", names(venn_intersections[i]), "** : ", file_name, "\n"))
}
message("Files with DEPs intersections were saved in \n", filenames)

Note:

Files with DEPs intersections were saved in

all : /home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/DEP_all.tsv

up : /home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/DEP_up.tsv

down : /home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/DEP_down.tsv

DEP expression divergences

We will remove from every intersection of rbohC Venn diagram the protein IDs that are UP or DOWN in the same intersection region in the other UP and DOWN Venn diagrams.

# define empty variables
proteins_divergent <- list()
discrep_prot <- data.frame(matrix(ncol = 2, nrow = 0))

# fill the variables with $all data
for (i in 1:length(venn_intersections[[i]])) {
  # remove from the intersection in $all those that are all UP DEPs
  tmp <- setdiff(venn_intersections$all[[i]], venn_intersections$up[[i]])
  # remove from remaining ids in the intersection in $all those that are all DOWN DEPs
  tmp <- setdiff(tmp, venn_intersections$down[[i]])
  proteins_divergent[[names(venn_intersections$all[i])]] <- tmp  # will be used below
  discrep_prot[nrow(discrep_prot) +1,] <- c(names(venn_intersections$all[i]),
                                            length(tmp))
}

The summary of divergences must be read as follows:

  • Intersection is the intersection region from the rbohC Venn diagram.
  • up refers to DEPs that are up-regulated in all genotypes compared in the intersection
  • down refers to DEPs that are down-regulated in all genotypes compared in the intersection
  • Divergent refers to DEPs that are down- or up-regulated in one genotype and up- or down-regulated in the others in the intersection.
# customize data to present as a kable
venn_data$intersections <- gsub(":", "-", venn_data$intersections) # change : to -
venn_data <- merge(discrep_prot, venn_data, by = 1, sort = FALSE)
colnames(venn_data)[1] <- "intersections"
colnames(venn_data)[2] <- "divergent"
venn_data <- venn_data[, -3]
venn_data <- venn_data[,c(1, 3, 4, 2)]

kable(venn_data,
      align = "lccc")
intersections up down divergent
rbohC-rbohF-WT 15 4 4
rbohC-rbohF 15 13 8
rbohC-WT 12 4 12
rbohF-WT 24 13 19
rbohC 121 138 0
rbohF 175 168 0
WT 216 180 0 
message("The normal situation is that no _divergent_ DEP is found in single comparisons")

Note:

The normal situation is that no divergent DEP is found in single comparisons

# remove useless variables
rm(tmp, tmp_up, tmp_down, name, intersections, proteins_ids, venn_data, discrep_prot)
rm(list=ls(pattern="^a[0-9]"))

Comparison of key DEPs

We represent a heatmap with the DEP that are common at least to two of the three conditions (WT, rbohC and rbohF) and are divergent.

proteins_intersections <- unlist(proteins_divergent, use.names = TRUE)

# prepare data for heatmap
heatmap_intersections <- list()
for (i in 1:length(counts)) {
  tmp <- counts[[i]][proteins_intersections,
                     grep("Abundance Ratio \\(log2\\)", colnames(counts[[i]])),
                     drop = FALSE]
  tmp <- na.omit(tmp)
  tmp$proteins <- rownames(tmp)
  tmp <- tmp[,c(2,1)]
  colnames(tmp)[2] <- "log2"
  heatmap_intersections[[names(counts[i])]] <- tmp
}

heatmap_intersections <- do.call(rbind, heatmap_intersections)
heatmap_intersections$condition <- rownames(heatmap_intersections)
heatmap_intersections$condition <- gsub("\\..*", "", heatmap_intersections$condition)

# represent heatmap
order <- as.data.frame(proteins_intersections)
order$intersections <- rownames(order)
order$intersections <- gsub("[0-9]", "", order$intersections)
order_name <- unique(order$intersections)

for (i in 1:length(order_name)) {
  order_proteins <- order$proteins_intersections[grep(paste0("^",
                                                       order_name[i],
                                                       "$"),
                                                order$intersections)]
  tmp2 <- data.frame()
  for (j in 1:length(order_proteins)) {
    tmp <- heatmap_intersections[grep(order_proteins[j],
                              heatmap_intersections$proteins),]
    tmp2 <- rbind(tmp2, tmp)
} # end for j
  heatmap_plot <- ggplot(tmp2,
                         aes(x = condition,
                             y = proteins,
                             fill = log2)) +
    geom_tile() + # rectangles with log2 values
    coord_equal() + # to make them squares
    scale_fill_gradient2(low = "#F8766D",
                         high = "#00BFC4",
                         mid = "white",
                         midpoint = 0,
                         na.value = "white") +
    scale_y_discrete(limits = rev) + # reverse order for y axis
    theme(#text = element_text(size = 5),
          # axis.text.x = element_text(angle = 90), # x labels vertical
          panel.grid.major = element_blank(), # remove major grid
          panel.grid.minor = element_blank(), # remove minor grid
          panel.background = element_blank()) + # remove background
    ggtitle(paste0("Divergent in ",
                   order_name[i]))
  print(heatmap_plot)
  } # end for i

Note that a white color indicates that the protein was not DEP for this condition.

# remove useless variables
rm(heatmap_intersections, heatmap_plot, tmp, proteins_intersections, order)

Common DEPs in rbohC and rbohF

We represent a heatmap with the DEP that are common to the two mutants (rbohC and rbohF) and are up or down.

proteins_common <- venn_intersections[c("up", "down")]
proteins_common_up <- proteins_common$up$`rbohC-rbohF`
proteins_common_down <- proteins_common$down$`rbohC-rbohF`
proteins_common <- c(proteins_common_up, proteins_common_down)

# prepare data for heatmap
heatmap_common <- list()
for (i in 1:length(counts)) {
  tmp <- counts[[i]][proteins_common,
                     grep("Abundance Ratio \\(log2\\)", colnames(counts[[i]])),
                     drop = FALSE]
  tmp <- na.omit(tmp)
  tmp$proteins <- rownames(tmp)
  tmp <- tmp[,c(2,1)]
  colnames(tmp)[2] <- "log2"
  heatmap_common[[names(counts[i])]] <- tmp
}

heatmap_common <- do.call(rbind, heatmap_common)
heatmap_common$condition <- rownames(heatmap_common)
heatmap_common$condition <- gsub("\\..*", "", heatmap_common$condition)

# represent heatmap
order <- heatmap_common[order(heatmap_common$log2,
                              decreasing = TRUE),]
order <- order[grep("rbohC", order$condition),]
order <- unique(order$proteins)
heatmap_common$proteins <- factor(heatmap_common$proteins,
                                  levels = order)

heatmap_plot <- ggplot(heatmap_common,
                       aes(x = condition,
                           y = proteins,
                           fill = log2)) +
  geom_tile() + # rectangles with log2 values
  coord_equal() + # to make them squares
  scale_fill_gradient2(low = "#F8766D",
                       high = "#00BFC4",
                       mid = "white",
                       midpoint = 0,
                       na.value = "white") +
  scale_y_discrete(limits = rev) + # reverse order for y axis
  theme(#text = element_text(size = 5),
        axis.text.x = element_text(angle = 90), # x labels vertical
        panel.grid.major = element_blank(), # remove major grid
        panel.grid.minor = element_blank(), # remove minor grid
        panel.background = element_blank()) # remove background
print(heatmap_plot)

# remove useless variables
rm(proteins_common, proteins_common_up, proteins_common_down, heatmap_common, heatmap_plot, tmp)

Transmembrane proteins

Let’s check how many of our DEPs are transmembrane.

transmembrane <- data.frame(matrix(ncol = 4, nrow = 0))
colnames(transmembrane) <- c("condition", "total", "trans", "perc")
for (i in 1:length(counts)) {
  tmp <- merge(x = counts[[i]],
               y = location,
               by = "row.names",
               all = FALSE,
               sort = FALSE)
  tmp <- tmp[c("Row.names", "Transmembrane")]
  total_prot <- length(rownames(tmp))
  keep <- grep("TRANSMEM", tmp$Transmembrane, ignore.case = TRUE)
  tmp <- tmp[keep,]
  transmembrane_prot <- length(rownames(tmp))
  transmembrane_perc <- round((transmembrane_prot/total_prot)*100,
                              digits = 0)
  
  transmembrane[nrow(transmembrane) +1,] <- c(names(counts[i]),
                                              total_prot,
                                              transmembrane_prot,
                                              transmembrane_perc)
  
  helical <- length(grep("Helical", tmp$Transmembrane, ignore.case = TRUE))
  not_helical <- FALSE
  if (helical != transmembrane_prot) {
    message(paste0("For ",
                 names(counts[i]),
                 " not all transmembrane proteins are helical type."))
    not_helical <- TRUE
  }
}

kable(transmembrane,
      col.names = c("Condition",
                    "Total",
                    "Transmemb.",
                    "Percent. (%)"),
      align = "lcccc")
Condition Total Transmemb. Percent. (%)
rbohC 322 80 25
rbohF 427 126 30
WT 468 133 28

Let’s see if there are transmembrane proteins that are not \(\alpha\)-helical.

if(not_helical == FALSE) {
  message("All transmembrane proteins are helical type.")
}

Note:

All transmembrane proteins are helical type.

# remove useless variables
rm(tmp, total_prot, transmembrane_prot, transmembrane_perc, helical, transmembrane)

Location distributions

Let’s see the distribution of DEPs by cellular location. Since there are many transmembrane proteins, two representations will be considered: one including all locations and another only for membrane locations. The last one will be analyzed with and without the non-specific membrane location.

All proteins, all locations

We show the all location for all proteins deteted in the experiment, being DEPs or not.

# assign location to all proteins
all_proteins_location <- merge(x = all_proteins,
                               y = location,
                               by = "row.names",
                               all = FALSE,
                               sort = FALSE)
all_proteins_location <- all_proteins_location[c("Row.names", "Subcellular.location")]
colnames(all_proteins_location)[1] <- "Protein.ID"

# split those proteins according to their location
subcell_all_proteins <- list()
data_barplot <- data.frame(matrix(ncol = 2, nrow = 0))
colnames(data_barplot) <- c("location", "count")
for (i in 1:length(location_cat)) {
  keep <- grep(location_cat[i],
               all_proteins_location$Subcellular.location,
               ignore.case = TRUE)
  tmp <- all_proteins_location[keep,]
  tmp <- tmp$Protein.ID
  empty <- length(tmp) == 0
  if (empty == FALSE) {
    subcell_all_proteins[[location_cat[i]]] <- tmp
    data_barplot[nrow(data_barplot) +1,] <- c(location_cat[i], length(tmp))
  }
}
data_barplot$count <- as.numeric(data_barplot$count)
# for membrane proteins
data_barplot_mem <- data_barplot[grep("membrane",
                                      data_barplot$location,
                                      ignore.case = TRUE),]
#
barplot_order <- data_barplot[order(data_barplot$count,
                                   decreasing = TRUE),]
barplot_order <- barplot_order$location
data_barplot$location <- factor(data_barplot$location,
                                levels = barplot_order)

# represent the results
barplot_plot <- ggplot(data_barplot,
                       aes(x = location,
                           y = count)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, # x label vertical 
                                   hjust = 1,
                                   vjust = 0.5)) +
  theme(panel.grid.major = element_blank(), # remove major grid
        panel.grid.minor = element_blank(), # remove minor grid
        panel.background = element_blank()) + # remove background
  ylim(0, 3000) # limits for y axis

ggplotly(barplot_plot)

It can be seen that the non-specific membrane location is the major location.

# remove useless variables
rm(tmp, all_proteins, all_proteins_location, subcell_all_proteins, data_barplot, barplot_order, barplot_plot)

All proteins, only membrane locations

barplot_order <- data_barplot_mem[order(data_barplot_mem$count,
                                        decreasing = TRUE),"location"]
data_barplot_mem$location <- factor(data_barplot_mem$location,
                                    levels = barplot_order)

# represent the results
barplot_plot <- ggplot(data_barplot_mem,
                       aes(x = location,
                           y = count)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, # x label vertical 
                                   hjust = 1,
                                   vjust = 0.5)) +
  theme(panel.grid.major = element_blank(), # remove major grid
        panel.grid.minor = element_blank(), # remove minor grid
        panel.background = element_blank()) + # remove background
  ylim(0, 3000) # limits for y axis

ggplotly(barplot_plot)
# remove useless variables
rm(barplot_order, data_barplot_mem, barplot_plot)

Membrane DEPs

We check the location of the membrane proteins of the comparisons made with the Venn diagram.

# assign location to the Venn diagram proteins
membrane_proteins_location <- list()
for (i in 1:length(venn_intersections)) {
  membrane_proteins_location[[names(venn_intersections[i])]] <- list()
  for (j in 1:length(venn_intersections[[i]])) {
    tmp <- venn_intersections[[i]][[j]]
    tmp <- merge(x = tmp,
                 y = location,
                 by.x = 1,
                 by.y = "row.names",
                 all = FALSE,
                 sort = FALSE)
    tmp <- tmp[c("x", "Subcellular.location")]
    colnames(tmp)[1] <- "Protein.ID"
    # keep only membrane proteins
    tmp <- tmp[grep("membrane",
                    tmp$Subcellular.location,
                    ignore.case = TRUE),]
    if(length(tmp) != 0) {
      name <- names(venn_intersections[[i]][j])
      membrane_proteins_location[[names(venn_intersections[i])]][[name]] <- tmp
    }
  }
}

# split those proteins according to their location
location_membrane <- location_cat[grep("membrane",
                                       location_cat,
                                       ignore.case = TRUE)]

subcell_membrane_proteins <- list()
subcell_save <- list()
for (i in 1:length(membrane_proteins_location)) {
  name <- names(membrane_proteins_location[i])
  subcell_membrane_proteins[[name]] <- data.frame(matrix(ncol = 3, nrow = 0))
  colnames(subcell_membrane_proteins[[name]]) <- c("location", "count", "intersections")
  subcell_save[[name]] <- data.frame(matrix(ncol = 3, nrow = 0))
  colnames(subcell_save[[name]]) <- c("protein", "intersections", "location")
  tmp <- membrane_proteins_location[[i]]
  tmp <- do.call(rbind, tmp)
  tmp$intersections <- rownames(tmp)
  tmp$intersections <- gsub("\\..*", "", tmp$intersections)
  for (j in 1:length(location_membrane)) {
    keep <- grep(location_membrane[j],
                 tmp$Subcellular.location,
                 ignore.case = TRUE)
    if (length(keep) != 0) {
      tmp2 <- tmp[keep,]
      tmp2_intersections <- tmp2$intersections
      tmp2_intersections <- unique(tmp2_intersections)
      for (k in 1:length(tmp2_intersections)) {
        tmp3 <- tmp2[grep(paste0("^", tmp2_intersections[k], "$"),
                          tmp2$intersections),]
        subcell_membrane_proteins[[i]][nrow(subcell_membrane_proteins[[i]]) +1,] <- c(location_membrane[j], length(rownames(tmp3)), tmp2_intersections[k])
        tmp3$location_membrane <- location_membrane[j]
        tmp3 <- tmp3[,-2]
        subcell_save[[i]] <- rbind(subcell_save[[i]], tmp3)
      } # end for k
    } # end if
  } # end for j
} # end for i

# for no membrane general location
subcell_membrane_no_proteins <- list()
for (i in 1:length(subcell_membrane_proteins)) {
  tmp <- subcell_membrane_proteins[[i]]
  tmp <- tmp[-grep("^membrane$", tmp$location, ignore.case = TRUE),]
  name <- names(subcell_membrane_proteins[i])
  subcell_membrane_no_proteins[[name]] <- tmp
}

# remove useless variables
rm(membrane_proteins_location, tmp, tmp2, tmp3, tmp2_intersections)

# with membrane general term
data_barplot <- list()
# represent the results
plt <- htmltools::tagList()
for (i in 1:length(subcell_membrane_proteins)) {
  data_barplot[[names(subcell_membrane_proteins[i])]] <- list()
  tmp <- subcell_membrane_proteins[[i]]
  tmp_intersections <- unique(tmp$intersections)
  l <- list()
  for (j in 1:length(tmp_intersections)) {
    l[[tmp_intersections[j]]] <- tmp[grep(paste0("^", tmp_intersections[j], "$"),
                                    tmp$intersections),c(1,2)]
  } # end for j
  df <- data.frame()
  for (k in 1:length(l)) {
    tmp2 <- l[[k]]
    colnames(tmp2)[2] <- tmp_intersections[k]
    if (k == 1){
      df <- tmp2
    } else {
      df <- merge(x = df,
                  y = tmp2,
                  by = 1,
                  all = TRUE,
                  sort = FALSE)
    } # end else
  } # end for k
  df[is.na(df)] <- 0
  mm <- melt(df, id.vars = "location")
  mm$value <- as.numeric(mm$value)
  
  gg_barplot<- ggplot(mm, aes(x = location, y = value)) +
    geom_bar(stat = "identity") + # can add width = 0.2 to change it
    facet_grid(.~variable) +
    coord_flip() +
    ggtitle(paste0(names(subcell_membrane_proteins[i]), " proteins")) +
    theme(panel.grid.major = element_blank(), # remove major grid
          panel.grid.minor = element_blank(), # remove minor grid
          panel.background = element_blank()) + # remove background
    theme(text = element_text(size = 8))
    
  plt[[i]] <- as_widget(ggplotly(gg_barplot))
} # end for i

Now, we represent the same results without the non-specific membrane location to highligth the abundance of DEPs in other membrane locations.

# no membrane general term
data_barplot <- list()
# represent the results
plt <- htmltools::tagList()
for (i in 1:length(subcell_membrane_no_proteins)) {
  data_barplot[[names(subcell_membrane_no_proteins[i])]] <- list()
  tmp <- subcell_membrane_no_proteins[[i]]
  tmp_intersections <- unique(tmp$intersections)
  l <- list()
  for (j in 1:length(tmp_intersections)) {
    l[[tmp_intersections[j]]] <- tmp[grep(paste0("^", tmp_intersections[j], "$"),
                                    tmp$intersections),c(1,2)]
  } # end for j
  df <- data.frame()
  for (k in 1:length(l)) {
    tmp2 <- l[[k]]
    colnames(tmp2)[2] <- tmp_intersections[k]
    if (k == 1){
      df <- tmp2
    } else {
      df <- merge(x = df,
                  y = tmp2,
                  by = 1,
                  all = TRUE,
                  sort = FALSE)
    } # end else
  } # end for k
  df[is.na(df)] <- 0
  mm <- melt(df, id.vars = "location")
  mm$value <- as.numeric(mm$value)
  
  gg_barplot<- ggplot(mm, aes(x = location, y = value)) +
    geom_bar(stat = "identity") + # can add width = 0.2 to change it
    facet_grid(.~variable) +
    coord_flip() +
    ggtitle(paste0(names(subcell_membrane_no_proteins[i]), " proteins")) +
    theme(panel.grid.major = element_blank(), # remove major grid
          panel.grid.minor = element_blank(), # remove minor grid
          panel.background = element_blank()) + # remove background
    theme(text = element_text(size = 8))
    
  plt[[i]] <- as_widget(ggplotly(gg_barplot))
} # end for i

Saving membrane DEPs.

Note:

Files with membrane DEPs were saved in

all : /home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/DEP_all_membrane.tsv

up : /home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/DEP_up_membrane.tsv

down : /home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/DEP_down_membrane.tsv

Summary table

A summary table with the protein ID, its subcellular location, if they are transmembrane (TRUE) or not (empty cell), the log2 value for the three conditions (rbohC, rbohF and WT) and the protein name.

datatable(final_datatable[, c(7, 1, 2, 3, 4, 5, 6)],  # reorder table
          options = list(pageLength = 5),  # number of lines by default
          caption = "Summary table of all the DEPs")
# save final_datatable
file_name <- paste0(SUB_DIR,
                    "summary_table.tsv")
write.table(final_datatable,
            file = file_name,
            quote = FALSE,
            sep = "\t",
            row.names = FALSE,
            col.names = TRUE)
message("The summary table was saved at \n", file_name, "\n")

Note:

The summary table was saved at

/home/bullones/Dropbox/Proteomica_arabidopsis/supplementary/results_2026-01-20_12.29.01/summary_table.tsv

About this session

Elapsed time: 0.06 min.

## Variables in memory:
##  [1] "col_min_max"                  "cols"                         "counts"                       "data_barplot"                
##  [5] "df"                           "distance"                     "empty"                        "FILE"                        
##  [9] "file_name"                    "FILE_NAME"                    "FILE_TOTAL"                   "filenames"                   
## [13] "final_datatable"              "gg_barplot"                   "heatmap_list"                 "HOY"                         
## [17] "i"                            "j"                            "k"                            "keep"                        
## [21] "l"                            "location"                     "LOCATION"                     "location_cat"                
## [25] "LOCATION_CAT"                 "location_membrane"            "mm"                           "my_palette"                  
## [29] "MY_TITLE"                     "name"                         "not_helical"                  "order"                       
## [33] "order_name"                   "order_proteins"               "plt"                          "proteins_divergent"          
## [37] "SOURCE_DIR"                   "SUB_DIR"                      "subcell_membrane_no_proteins" "subcell_membrane_proteins"   
## [41] "subcell_save"                 "T_total"                      "T00"                          "tmp"                         
## [45] "tmp_hc"                       "tmp_intersections"            "tmp_venn"                     "tmp2"                        
## [49] "venn_intersections"           "venn_plot"                    "VERBOSE_MODE"                 "WD_DIR"                      
## 
## R version 4.5.2 (2025-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=es_ES.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=es_ES.UTF-8     LC_MONETARY=es_ES.UTF-8   
##  [6] LC_MESSAGES=es_ES.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Europe/Madrid
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] DT_0.34.0           reshape2_1.4.5      data.table_1.17.8   dplyr_1.1.4         plotly_4.11.0       VennDiagram_1.7.3   futile.logger_1.4.3
##  [8] ggplot2_4.0.1       venn_1.12           rmarkdown_2.30      knitr_1.50         
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10          generics_0.1.4       tidyr_1.3.1          futile.options_1.0.1 stringi_1.8.7        digest_0.6.38       
##  [7] magrittr_2.0.4       evaluate_1.0.5       RColorBrewer_1.1-3   fastmap_1.2.0        plyr_1.8.9           jsonlite_2.0.0      
## [13] formatR_1.14         httr_1.4.7           purrr_1.2.0          crosstalk_1.2.2      viridisLite_0.4.2    scales_1.4.0        
## [19] lazyeval_0.2.2       jquerylib_0.1.4      cli_3.6.5            rlang_1.1.6          cachem_1.1.0         withr_3.0.2         
## [25] yaml_2.3.10          tools_4.5.2          lambda.r_1.2.4       vctrs_0.6.5          R6_2.6.1             lifecycle_1.0.4     
## [31] stringr_1.6.0        htmlwidgets_1.6.4    admisc_0.39          pkgconfig_2.0.3      pillar_1.11.1        bslib_0.9.0         
## [37] gtable_0.3.6         rsconnect_1.6.1      glue_1.8.0           Rcpp_1.1.0           xfun_0.54            tibble_3.3.0        
## [43] tidyselect_1.2.1     rstudioapi_0.17.1    dichromat_2.0-0.1    farver_2.1.2         htmltools_0.5.8.1    labeling_0.4.3      
## [49] compiler_4.5.2       S7_0.2.1